#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>
#include <gmp.h>
#include <pthread.h>

#define BIT_PRECISION 33219281

struct thread_data {
	mpf_t a0;
	mpf_t a1;
	mpf_t y0;
	mpf_t y1;
	mpf_t pi;
	int iteration;
};

pthread_t thread_y, thread_a;
pthread_mutex_t mutex;

mpf_t prec;
mpf_t sqroot;
int conv;

FILE *fp;

void *calculate_y(void *args) {
	pthread_mutex_lock(&mutex);

	mpf_t aux, aux2, aux3;

	struct thread_data *data = args;

	mpf_init(aux);
	mpf_init(aux2);
	mpf_init(aux3);

	// y0 = y0 ^ 4
	mpf_pow_ui(data->y0, data->y0, 4);
	// aux = 1 - y0
	mpf_ui_sub(aux, 1, data->y0);
	// aux = Raiz(aux)
	mpf_sqrt(aux, aux);
	// aux = Raiz(aux);
	mpf_sqrt(aux, aux);
	// aux2 = 1 - aux
	mpf_ui_sub(aux2, 1, aux);
	// aux3 = 1 + aux
	mpf_add_ui(aux3, aux, 1);
	// y1 = aux2 / aux3
	mpf_div(data->y1, aux2, aux3);

	mpf_clear(aux);
	mpf_clear(aux2);
	mpf_clear(aux3);

	pthread_mutex_unlock(&mutex);
	pthread_exit(NULL);
}

void *calculate_a(void *args) {
	pthread_mutex_lock(&mutex);

	mpf_t aux, aux2, aux3;
	unsigned long int power_of_two;

	struct thread_data *data = args;

	mpf_init(aux);
	mpf_init(aux2);
	mpf_init(aux3);

	// aux = y1 ^ 2
	mpf_pow_ui(aux, data->y1, 2);
	// aux = aux + y1
	mpf_add(aux, aux, data->y1);
	// aux = aux + 1
	mpf_add_ui(aux, aux, 1);

	// power_of_two = 2 * i + 3
	power_of_two = pow(2, (2 * data->iteration) + 3);
	// aux2 = y1 * power_of_two
	mpf_mul_ui(aux2, data->y1, power_of_two);

	// aux3 = aux * aux2
	mpf_mul(aux3, aux, aux2);
	//mpf_out_str(stdout, 10, 0, aux2);

	// aux = y1 + 1
	mpf_add_ui(aux, data->y1, 1);
	// aux = aux ^ 4
	mpf_pow_ui(aux, aux, 4);
	// aux = aux * a0
	mpf_mul(aux, aux, data->a0);

	// a1 = aux - aux3
	mpf_sub(data->a1, aux, aux3);

	// Verifica a precisao
	mpf_sub(aux, data->a1, data->a0);
	mpf_abs(aux, aux);
	conv = mpf_cmp(aux, prec);

	mpf_set(data->a0, data->a1);
	mpf_set(data->y0, data->y1);
	mpf_ui_div(data->pi, 1, data->a1);
	
	mpf_clear(aux);
	mpf_clear(aux2);
	mpf_clear(aux3);

	pthread_mutex_unlock(&mutex);
	pthread_exit(NULL);
}

int main(void)
{
	struct timeval start, end; // Variaveis para guardar o tempo de execucao
	int time;

	mpf_set_default_prec(BIT_PRECISION);

	mpf_t a0, a1, y0, y1, prec;
	mpf_t aux, aux2, aux3;
	struct thread_data data;

	fp = fopen("borwein_paralelo.out", "w+");

	pthread_mutex_init(&mutex, NULL);

	mpf_init(data.a0);
	mpf_init(data.a1);
	mpf_init(data.y0);
	mpf_init(data.y1);
	mpf_init(data.pi);

	mpf_init_set_ui(a0, 6);
	mpf_init(a1);
	mpf_init(y0);
	mpf_init(y1);
	mpf_init(aux);
	mpf_init_set_ui(aux2, 1);
	mpf_init(aux3);
	mpf_init(sqroot);
	mpf_init(prec);

	printf("Calculando precisao...\n");

	// Calculo da precisao (condicao de parada)
	// prec = 1 / 10 ^ 10 000 000
	mpf_init_set_ui(aux2, 1);
	mpf_init_set_ui(prec, 10);
	mpf_pow_ui(prec, prec, 10000000);
	mpf_ui_div(prec, 1, prec);

	gettimeofday(&start, NULL);
	printf("Calculando pi...\n");

	// Raiz(2)
	mpf_sqrt_ui(sqroot, 2);
	// aux = 4 * Raiz(2)
	mpf_mul_ui(aux, sqroot, 4);
	// a0 = 6 - 4 * Raiz(2)
	mpf_sub(a0, a0, aux);

	mpf_set(data.a0, a0);

	// y0 = Raiz(2) - 1
	mpf_sub_ui(y0, sqroot, 1);

	mpf_set(data.y0, y0);

	data.iteration = 0;

	int i;
	for (i = 0; conv > 0 || i == 0; i++) {
		pthread_create(&thread_y, NULL, &calculate_y, &data);
		pthread_join(thread_y, NULL);
		pthread_create(&thread_a, NULL, &calculate_a, &data);
		pthread_join(thread_a, NULL);
		data.iteration++;
	}

	printf("Pronto.\n");

	gettimeofday(&end, NULL);
	time = end.tv_sec - start.tv_sec;

	printf("\nNumero de iteracoes: %d\nTempo para calculo: %d segundos\n", data.iteration, time);

	printf("Escrevendo saida no arquivo...\n");
	mpf_ui_div(data.a0, 1, data.a0);
	mpf_out_str(fp, 10, 0, a0);
	printf("Pronto: borwein_paralelo.out\n");

	mpf_clear(data.a0);
	mpf_clear(data.a1);
	mpf_clear(data.y0);
	mpf_clear(data.y1);
	mpf_clear(a0);
	mpf_clear(a1);
	mpf_clear(y0);
	mpf_clear(y1);
	mpf_clear(aux);
	mpf_clear(aux2);
	mpf_clear(aux3);
	mpf_clear(sqroot);
	mpf_clear(prec);

	return EXIT_SUCCESS;
}
